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Abstract 

Computational methods are proposed for solving a convex quadratic program 
(QP). Active-set methods are defined for a particular primal and dual formu¬ 
lation of a QP with general equality constraints and simple lower bounds on 
the variables. In the first part of the paper, two methods are proposed, one 
primal and one dual. These methods generate a sequence of iterates that are 
feasible with respect to the equality constraints associated with the optimality 
conditions of the primal-dual form. The primal method maintains feasibility 
of the primal inequalities while driving the infeasibilities of the dual inequal¬ 
ities to zero. The dual method maintains feasibility of the dual inequalities 
while moving to satisfy the primal inequalities. In each of these methods, the 
search directions satisfy a KKT system of equations formed from Hessian and 
constraint components associated with an appropriate column basis. The com¬ 
position of the basis is specified by an active-set strategy that guarantees the 
nonsingularity of each set of KKT equations. Each of the proposed methods is 
a conventional active-set method in the sense that an initial primal- or dual- 
feasible point is required. In the second part of the paper, it is shown how 
the quadratic program may be solved as a coupled pair of primal and dual 
quadratic programs created from the original by simultaneously shifting the 
simple-bound constraints and adding a penalty term to the objective function. 

Any conventional column basis may be made optimal for such a primal-dual 
pair of shifted-penalized problems. The shifts are then updated using the solu¬ 
tion of either the primal or the dual shifted problem. An obvious application 
of this approach is to solve a shifted dual QP to define an initial feasible point 
for the primal (or vice versa). The computational performance of each of the 
proposed methods is evaluated on a set of convex problems from the CUTEst 
test collection. 
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Primal and Dual Methods for Convex QP 


1. Introduction 

We consider the formulation and analysis of active-set methods for a convex quadratic 
program (QP) of the form 


minimize 


Ix'^Hx + \y^My + c^x 


snbject to Ax + My = b, x > 0, 


( 1 . 1 ) 


where A, b, c, H and M are constant, with H and M symmetric positive semidefi- 


nite. In order to simplify the theoretical discussion, the inequalities of (1.1) involve 
nonnegativity constraints only. However, the methods to be described are easily 
extended to treat all forms of linear constraints. (Numerical results are given for 
problems with constraints in the form x^ < x < Xu and b^ < Ax < bu, for fixed vec¬ 
tors Xl, Xu, bu and bu-) If M = 0, the QP (1.1) is a conventional convex quadratic 


program with constraints defined in standard form. A regularized quadratic pro¬ 
gram may be obtained by defining M = fil for some small positive parameter y. 
(For applications that require the solution of a regularized QP see, e.g., (1 g[^.) 
Active-set methods for quadratic programming problems of the form (1.1) solve 


a sequence of linear equations that involve the y-variables and a snbset of the x- 
variables. Each set of equations constitutes the optimality conditions associated with 
an equality-constrained quadratic subproblem. The goal is to predict the optimal 
active set, i.e., the set of constraints that are satisfied with equality, at the solution 
of the problem. A conventional active-set method has two phases. In the hrst phase, 
a feasible point is found while ignoring the objective function; in the second phase, 
the objective is minimized while feasibility is maintained. A useful feature of active- 
set methods is that they are well-suited for “warm starts”, where a good estimate 
of the optimal active set is used to start the algorithm. This is particularly useful in 
applications where a sequence of quadratic programs is solved, e.g., in a sequential 
quadratic programming method or in an ODE- or PDE-constrained problem with 
mesh refinement. Other applications of active-set methods for quadratic program¬ 
ming include mixed-integer nonlinear programming, portfolio analysis, structural 
analysis, and optimal control. 

In Section the primal and dual forms of a convex quadratic program with 
constraints in standard form are generalized to include general lower bounds on both 
the primal and dual variables. These problems constitute a primal-dual pair that 
includes problem ( |1.1[ ) and its associated dual as a special case. In Sections and 
an active-set method is proposed for each of the primal and dual forms associated 
with the generalized problem of Section Both of these methods provide a sequence 
of iterates that are feasible with respect to the equality constraints associated with 
the optimality conditions of the primal-dual problem pair. The primal method 
maintains feasibility of the primal inequalities while driving the infeasibilities of the 
dual inequalities to zero. By contrast, the dual method maintains feasibility of the 
dual inequalities while moving to satisfy the primal inequalities. In each of these 
methods, the search directions satisfy a KKT system of equations formed from 
Hessian and constraint components associated with an appropriate column basis. 
The composition of the basis is specified by an active-set strategy that guarantees 
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the nonsingularity of each set of KKT equations. 

The methods formulated in Sections [SHU define conventional active-set methods 
in the sense that an initial feasible point is required. In Section a method is 
proposed that solves a pair of coupled quadratic programs created from the original 
by simultaneously shifting the simple-bound constraints and adding a penalty term 
to the objective function. Any conventional column basis can be made optimal for 
such a primal-dual pair of shifted-penalized problems. The shifts are then updated 
using the solution of either the primal or the dual shifted problem. An obvious 
application of this idea is to solve a shifted dual QP to define an initial feasible 
point for the primal, or vice-versa. In addition to the obvious benefit of using the 
objective function while getting feasible, this approach provides an effective method 
for finding a dual-feasible point when H is positive semidefinite and M = 0. Finding 
a dual-feasible point is relatively straightforward for the strictly convex case, i.e., 
when H is positive definite. However, in the general case, the dual constraints for 
the phase-one linear program involve entries from H as well as A, which complicates 
the formulation of the phase-one method considerably. 

Finally, in Section some numerical experiments are presented for a simple 
Matlab implementation of a coupled primal-dual method applied to a set of convex 


problems from the CUTEst test collection 143, 45 


There are a number of alternative active-set methods available for solving a QP 
with constraints written in the format of problem (1.1). Broadly speaking, these 
methods fall into three classes defined here in the order of increasing generality: 
(i) methods for strictly convex quadratic programming (H symmetric positive defi¬ 
nite) [^ [M|[4H(55|[58] ; (ii) methods for convex quadratic programming {H symmetric 
positive semidefinite) [8,35,^52,^; and (hi) methods for general quadratic pro¬ 
gramming (no assumptions on H other than symmetry) p|[4| [IH[^[^[^[37t[Mt 
|^|^|^|^[^[^. Of the methods specifically designed for convex quadratic 
programming, only the methods of Boland and Wong [59[ Chapter 4] are dual 
active-set methods. Some existing active-set quadratic programming solvers include 


QPOPT [33], qPSchur (^, SQOPT (^, SQIC (4^ and QPA (part of the GALAHAD software 
library)^ |44] . 

The primal active-set method proposed in Section is motivated by the methods 
of Fletcher [24] , Gould , and Gill and Wong , which may be viewed as methods 
that extend the properties of the simplex method to general quadratic programming. 
At each iteration, a direction is computed that satisfies a nonsingular system of 
linear equations based on an estimate of the active set at a solution. The equations 
may be written in symmetric form and involve both the primal and dual variables. 
In this context, the purpose of the active-set strategy is not only to obtain a good 
estimate of the optimal active set, but also to ensure that the systems of linear 
equations that must be solved at each iteration are nonsingular. This strategy 
allows the application of any convenient linear solver for the computation of the 
iterates. In this paper, these ideas are applied to convex quadratic programming. 
The resulting sequence of iterates is the same as that generated by an algorithm for 
general QP, but the structure of the iteration is different, as is the structure of the 
linear equations that must be solved. Similar ideas are used to formulate the new 
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dual active-set method proposed in Section 

The proposed primal, dual, and combined primal-dual methods use a “con¬ 
ventional” active-set approach in the sense that the constraints remain unchanged 
during the solution of a given QP. Alternative approaches that use a parametric 
active-set method have been proposed by Best [^|^, Ritter 56,57 , Ferreau, Bock 
and Diehl 1^, Potschka et al. [54], and implemented in the qpOASES package by Fer¬ 


reau et al. [^. Primal methods based on the augmented Lagrangian method have 
been proposed by Delbos and Gilbert [^, Chiche and Gilbert [^, and Gilbert 
and Joannopoulos 30 . The use of shifts for the bounds have been suggested by 


Cartis and Gould [13] in the context of interior methods for linear programming. 
Another class of active-set methods that are convergent for strictly convex quadratic 
programs have been considered by Gurtis, Han, and Robinson (TBl. 


Notation and terminology. Given vectors a and b with the same dimension, 
min(a, 6) is a vector with components min(ai,6j). The vectors e and ej denote, 
respectively, the column vector of ones and the jth column of the identity matrix 
I. The dimensions of e, e* and I are defined by the context. Given vectors x and y, 
the column vector consisting of the components of x augmented by the components 
of y is denoted by {x,y). 


2. Background 


Although the purpose of this paper is the solution of quadratic programs of the form 
(1.1), for reasons that will become evident in Section]^ the analysis will focus on the 
properties of a pair of problems that may be interpreted as a primal-dual pair of QPs 


associated with problem (1.1). It is assumed throughout that the matrix (A M) 
associated with the equality constraints of problem (1.1) has full row rank. This 


assumption can be made without loss of generality, as shown in Proposition A.8 of 
the Appendix. The paper involves a number of other basic theoretical results that 
are subsidiary to the main presentation. The proofs of these results are relegated 
to the Appendix. 


2.1. Formulation of the primal and dual problems 

For given constant vectors q and r, consider the pair of convex quadratic programs 

hx^Hx + hy'^My + c^x + r'^x 


and 


(PQPq,r) 

minimize 

x,y 

subject to 

(DQPg,Q 

maximize 

x,y,z 

subject to 


Ax + My = b, 


X > —q, 


— \x'^Hx — ^y'^My -|- b'^y — q'^z 

— Hx + A^y + z = c, z > —r. 


The following result gives joint optimality conditions for the triple (x, y, z) such 
that (x, y) is optimal for (PQPg ,,), and {x,y,z) is optimal for (DQPg ,,). If q and r 
are zero, then (PQPo,o) (DQPq g) are the primal and dual problems associated 
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with (1.1). For arbitrary q and r, (PQPg,,) and (DQP^,.) are essentially the dual of 
each other, the difference is only an additive constant in the value of the objective 
function. 

Proposition 2.1. Let q and r denote constant vectors in ML. If (x, y, z) is a given 
triple in M"’ x x then {x, y) is optimal for {PQPq .f) and {x, y, z) is optimal 
for (DQPq j.) if and only if 


o' 

II 

1 

1 

o 

+ 

so 

(2.1a) 

Ax + My — b = 0, 

(2.1b) 

x + q > 0, 

(2.1c) 

z + r > 0, 

(2.1d) 

(x + q)'^(z + r) = 0. 

(2.1e) 


In addition, the opt imal objective values satisfy optval(PQPq,.) — optval{DQPg j.) = 
—q^r. Finally, (2.1) has a solution if and only if the sets 


{(x, y, z) : —Hx + A^y + z = c, z > —r} and {(x, y) : Ax + My = b, x > —o'} 
are both nonempty. 

Proof. Let the vector of Lagrange multipliers for the constraints Ax + My — 6 = 0 
be denoted by y. Without loss of generality, the Lagrange multipliers for the bounds 
X + y > 0 of (PQPq,.) may be written in the form z + r, where r is the given fixed 
vector r. With these definitions, a Lagrangian function L(x, y, y, z) associated with 
(PQPg,r) is given by 

L(x, y, y, z) = \x^Hx + (c + r^x + \y'^My - y^{Ax + My - b) 

- (z + r)'^{x + q). 

Stationarity of the Lagrangian with respect to x and y implies that 

Hx + c + r — A^y — z — r = Hx + c — A^y — z = 0, (2.2a) 

My-My = 0. (2.2b) 

The optimality conditions for (PQPg r) then given by: (i) the feasibility condi¬ 


tions (2.1b) and (2.1c); (ii) the nonnegativity conditions (2.Id) for the multipliers 


associated with the bounds x + q > 0] (hi) the stationarity conditions (2.2); and (iv 


the complementarity conditions (2.1e). The vector y appears only in the term My 


of (2.1b) and (2.2b). In addition, (2.2b) implies that My = My, in which case we 


may choose y = y. This common value of y and y must satisfy (2.2a), which is then 


equivalent to (2.1a). The optimality conditions (2.1) for (PQP^,,) follow directly. 

With the substitution y = y, the expression for the Lagrangian may be rear¬ 
ranged so that 


L(x, y, y, z) = — \x^Hx — ^y'^My + b^y — q^z -|- {Hx + c — A^y — zY'x — q^r. (2.3) 
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Taking into account (2.2) for y = y, the dual objective is given by (2.3) as — \x'^Hx — 
^y'^My + h'^y — q^z — q^r, and the dual constraints are Hx + c — A^y — z = 0 and 
z + r > 0. It follows that (DQP^,.) is equivalent to the dual of (PQPg,.), the only 
difference is the constant term —q^r in the objective, which is a consequence of the 
shift 2 + r in the dual variables. Consequently, strong duality for convex quadratic 
programming implies optval(PQP^) — optval(DQPq,,) = —In addition, the 

and (DQPg ,,) with 
is 


variables x, y and z satisfying (2.1) are feasible for (PQPg 
the difference in the objective function value being —q^r. It follows that (x, y, z) 




optimal for (DQP^ ,.) as well as (PQP^ ,.). Finally, feasibility of both (PQPg ,,) and 
(DQPg is both necessary and sufficient for the existence of optimal solutions. | 


2.2. Optimality conditions and the KKT equations 

The proposed methods are based on maintaining index sets B and J\f that define 
a partition of the index set X = {1, 2, ..., n}, i.e., Z = B Li Af with B C M = %. 
Following standard terminology, we refer to the subvectors Xg and x^ associated 
with an arbitrary x as the basic and nonbasic variables, respectively. The crucial 
feature of B is that it dehnes a unique solution (x, y, z) to the equations 

Hx + c - A^y - z = Xiv + giv = 0, 

Ax + My — 6 = 0, Zg + Tg = Q. 


For the symmetric Hessian H, the matrices Hgg and denote the subset of rows 
and columns of H associated with the sets B and AA, respectively. The unsymmetric 
matrix of components hij with i £ B and j G Af will be denoted by Hgj^. Similarly, 
Ag and Aj^ denote the matrices of columns of A associated with B and Af respec¬ 


tively. With this notation, the equations (2.4) may be written in partitioned form 


as 


a^iv + g'jv = 0, 
Zg+rg = 0, 


HggXg + HgffXj^ + Cg - A^y - Zg=Q, 

^BN^B T Hj^j^X^ T Cjv -^nV ^5 

AgXg + AjvXjv + My -6 = 0. 

Eliminating Xjv and Zg from these equations using the equalities Xn + q^ = 0 and 
Zg + Vg = 0 yields the symmetric equations 


Hgg Al 
Ag -M 


Xg 

-y 


Hgf/qf^ — Cg — rg 

A^qN + 6 


(2.5) 


for Xg and y. It follows that (2.4) has a unique solution if and only if (2.5) has a 


unique solution. Therefore, if B is chosen to ensure that (2.4) has a unique solution. 


it must follow from (2.5) that the matrix Kg such that 


Kg = 


H 


Al 


BB 

-M 


( 2 . 6 ) 


is nonsingular. Once Xg and y have been computed, the z^v-variables are given by 

Zm Bgj^Xg — Hj^j^qj^ Cff — Aj^y. (2-7) 
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As in Gill and Wong |40] , any set B such that is nonsingular is referred to as 
a second-order eonsistent basis. Methods that impose restrictions on the eigenval¬ 
ues of Kb are known as inertia-controlling methods. (For a description of inertia¬ 
controlling methods for general quadratic programming, see, e.g.. Gill et al. [38] , 
and Gill and Wong (40].) 

The two methods proposed in this paper, one primal, one dual, generate a se¬ 


quence of iterates that satisfy the equations (2.4) for some partition B and Af. If the 
conditions ( |2.4[ ) are satisfied, the additional requirement for fulfilling the optimality 
conditions of Proposition [tT] are Xb+Qb >0 and > 0. The primal method of 

Section imposes the restriction that Xb + Qb > 0, which implies that the sequence 
of iterates is primal feasible. In this case the method terminates when > 0 

is satisfied. Conversely, the dual method of Section imposes dual feasibility by 
means of the bounds Zn + r^f > 0 and terminates when Xb + Qb > 0. 

In both methods, an iteration starts and ends with a second-order consistent 
basis, and comprises one or more subiterations. In each subiteration an index I and 
index sets B and Af are known such that .BU {/} UAf = {1, 2, ..., n}. This partition 
defines a search direction (Ax, Ay, Az) that satishes the identities 


HAx — A^Ay — Az = 0, Axn = 0, 

AAx + MAy = 0, Azb = 0. 


( 2 . 8 ) 


As I ^ B and I 0 Af, these conditions imply that neither Axi nor 
to be zero. The conditions Ax^ = 0 and Azb = 0 imply that (2.8) 
in the partitioned-matrix form 


Azi are restricted 
may be expressed 


fhii 

hi 


hsi 

HbB 

Af 

hm 

hL 

-^N 

\ ai 

Ab 

-M 



/ Axi\ 
Axb 
-A y 
-Azi 

\ — AZ]sr) 


/o\ 

0 

0 

Vo/ 


where hu denotes the /th diagonal of H, and the column vectors hgi and /ijv! denote 
the column vectors of elements hu and hji with i £ B, and j G Af, respectively. It 
follows that Axi, Axb, Ay and Azi satisfy the homogeneous equations 


(K hi af 1 

I hgi HBB -^b 

\ ai Ab -M 


( Axi\ 
Axb 
-A y 

\ / 



and Azjq is given by 


Azf^ hffiAxi -|- Hg^Ax B Aj^Ay. 


(2.9a) 


(2.9b) 


The properties of these equations are established in the next subsection. 
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2.3. The linear algebra framework 

This section establishes the linear algebra framework that serves to emphasize the 
underlying symmetry between the primal and dual methods. It is shown that the 
search direction for the primal and the dual method is a nonzero solution of the 


homogeneous equations (2.9a), i.e., every direction is a nontrivial null vector of the 


matrix of (2.9a). In particular, it is shown that the null-space of (2.9a) has dimension 


one, which implies that the solution of (2.9a) is unique up to a scalar multiple. The 


length of the direction is then completely determined by hxing either Axi = 1 or 
Azi = 1. The choice of which component to hx depends on whether or not the 
corresponding component in a null vector of (2.9a) is nonzero. The conditions are 


stated precisely in Propositions 2.3 and 2.4 below. 

The hrst result shows that the components Axi and Azi of any direction {Ax, 


Ay, Az) satisfying the identities (2.8) must be such that AxiAzi > 0. 


Proposition 2.2. If the vector {Ax, Ay, Az) satisfies the identities 

HAx — A^Ay — Az = 0, 

AAx + MAy = 0, 

then Ax'^Az = Ax^HAx -p Ay'^MAy > 0. Moreover, given an index I and index 
sets B and M such that B U {/} U AC = {1, 2, ..., n} with Axn = 0 and Azb = 0, 
then AxiAzi = Ax'^HAx + Ay'^MAy > 0. 

Proof. Premultiplying the hrst identity by Ax'^ and the second by Ay'^ gives 

Ax^H Ax — Ax^A^Ay — Ax'^Az = 0, and Ay'^AAx + Ay^M Ay = 0. 

Eliminating the term Ax'^A'^Ay gives Ax'^HAx -p Ay'^MAy = Ax'^Az. By dehni- 
tion, H and M are symmetric positive semidehnite, which gives Ax'^Az > 0. In 
particular, if .6 U {/} U AC = {1, 2, ..., n}, with Ax^ = 0 and Azb = 0, it must hold 
that Ax'^Az = AxiAzi >0. | 


The set of vectors {Axi, Axb, Ay, Azi, Aziq) satisfying the equations (2.9) is 


completely characterized by the properties of the matrices Kb and AT; such that 


= 




A'^ \ 

(ha 

hi 

ar 

-M J 

and Ki = 

Hbb 


\ ai 


-M 


( 2 . 10 ) 


The properties are summarized by the results of the following two propositions. 

Proposition 2.3. Assume that Kb is nonsingular. Let Axi be a given nonnegative 
scalar. 


1. If Ax I = 0, then the only solution of (2.9) is zero, i.e., Axb = 0, Ay = 0, 
Azi = 0 and Az^ = 0. 
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2. If Ax I > 0, then the quantities Axb, Ay, Azi and Az^ of {2.9) are unique and 


satisfy the equations 


Ab -m 


^ _ fhgi 
-Ay ) \ai 

Azi = hiiAxi + hl^Axg - afAy, 

Azff -|- Hgj^Axg Aj^Ay. 


Axi, 

.T 

T 


( 2 . 11 ) 


Moreover, either 

(i) Ki is nonsingular and Azi > 0, or 

(ii) Ki is singular and Azi = 0, in which case it holds that Ay = 0, Az^ = 0, 
and the multiplicity of the zero eigenvalue of Ki is one, with corresponding 
eigenvector {Axi, Axb,0)- 


Proof. Proposition 2.2 implies that Azi > 0 if Axi > 0, which implies that the 
statement of the proposition includes all possible values of Azi. The second and 


third blocks of the equations (2.9a) imply that 


hgi 

ai 


Axi + 


Hbb 

Ab 


-M 


Axb 

-Ay 


( 2 . 12 ) 


As Kg is nonsingular by assumption, the vectors Axg and Ay must constitute the 


unique solution of (2.12) for a given value of Axi. Furthermore, given Axg and Ay, 


the quantities Azi and Az^ of (2.11) are also uniquely defined. The specific value 
Axi = 0, gives Axb = 0 and Ay = 0, so that Azi = 0 and Azn = 0. It follows 
that Axi must be nonzero for at least one of the vectors Axg, Ay, Azi or Az^ to 
be nonzero. 

Next it is shown that if Axi > 0, then either ( [^ or (2ii) must hold. For ([^, it 
is necessary to show that if Axi > 0 and Ki is nonsingular, then Azi > 0. If Ki is 


nonsingular, the homogeneous equations (2.9a) may be written in the form 




(2.13) 


which implies that Axi, Axg and Ay are unique for a given value of Azi. In 
particular, if Azi = 0 then Axi = 0, which would contradict the assumption that 
Axi > 0. If follows that Azi must be nonzero. Finally, Proposition |2.2| implies that 
if Azi is nonzero and Axi > 0, then Azi > 0 as required. 


For the first part of (2ii), it must be shown that if Ki is singular, then Azi = 0. 


If Ki is singular, it must have a nontrivial null vector (pi, pg, —u). Moreover, every 
null vector must have a nonzero pi, because otherwise {pg, —u) would be a nontrivial 
null vector oi Kg, which contradicts the assumption that Kg is nonsingular. A fixed 
value of Pi uniquely defines pg and u, which indicates that the multiplicity of the 
zero eigenvalue must be one. A simple substitution shows that {pi, pg, —u, vi) 
is a nontrivial solution of the homogeneous equation (2.9a) such that vi = 0. As 
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the subspace of vectors satisfying (2.9a) is of dimension one, it follows that every 


solution is unique up to a scalar multiple. Given the properties of the known solution 


{Ph Pb, —u, 0), it follows that every solution {Axi, Axb, —Ay, —Azi) of (2.9a) is 
an eigenvector associated with the zero eigenvalue of Ki, with Azi = 0. 


For the second part of (2ii), if Azi = 0, the homogeneous equations (2.9a) become 


h^ 

"'Bl 

of \ 

/ Axi 

Hgg 

Al 

AXg 

Ag 

-Mj 

\-Ay 


(2.14) 


As Ki is singular in (2.14), Proposition A.l of the Appendix implies that 
'h,. 



and 



Ay = 


(2.15) 


The nonsingularity of Kg implies that (Ag — M) has full row rank, in which case 
the second equation of (2.15) gives Ay = 0. It follows that every eigenvector of AT; 


associated with the zero eigenvalue has the form {Axi, Axg, 0). It remains to show 
that Azn = 0. If Proposition |A.2| of the Appendix is applied to the first equation 


of (2.15), then it must hold that 


Hg 

Hi 


( Axi 

\AXg 


It follows from the definition of Azj^ in (2.11) that Azj^ = h^^^Axi + H^j^Axg — 
A^^Ay = 0, which completes the proof. | 


Proposition 2.4. Assume that Ki is nonsingular. Let Azi be a given nonnegative 
scalar. 


1. If Azi = 0, then the only solution of (2.9) is zero, i.e., Axi = 0, Axg = 0, 
Ay = 0 and Az^ = 0. 


2. If Azi > 0, then the quantities Axi, Axg, Ay and Az^ of {2.9) are unique and 
satisfy the equations 

Axg I = 1^0 j Azi, (2.16a) 

Azff = Hj^i^Axi + Hg^Axg - A^Ay. (2.16b) 


hm 

\ 

( Axi' 

Hgg 

Al 

AXg 

Ag 

-MJ 

\-Ay , 


Moreover, either 

(i) Kg is nonsingular and Axi > 0, or 

(ii) Kg is singular and Axi = 0, in which case, it holds that Axg = 0 and 
the multiplicity of the zero eigenvalue of Kg is one, with corresponding 
eigenvector (0, Ay). 
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Proof. In Proposition 2.2 it is established that Axi > 0 if Azi > 0, which implies 
that the statement of the proposition includes all possible values of Axi. 

It follows from (2.9a) that Axi, Axg, and Ay must satisfy the equations 


hi 

\ 

/ Axi 

Hbb 


Axb 

Ab 

-mJ 

\-^y 


(2.17) 


Under the given assumption that Ki is nonsingular, the vectors Axi, Axb and Ay 



are uniquely determined by (2.17) for a hxed value of Azi. In addition, once Axi^ 
Axb and Ay are dehned, Azj^ is uniquely determined by (2.16b). It follows that if 
Azi = 0, then Axi = 0, Axb = 0, Ay = 0 and Az^^ = 0. 

It remains to show that if Azi > 0, then either 


i^ or (2ii) must hold. If Kb is 


singular, then Proposition |A.l of the Appendix implies that there must exist u and 
V such that 


Proposition A.2 of the Appendix implies that the vector u must also satisfy hgi'^ — 0- 


If u is nonzero, then (0, u, 0) is a nontrivial null vector for Ki, which contradicts the 
assumption that Ki is nonsingular. It follows that ( Hbb Ag ) has full row rank 
and the singularity of Kb must be caused by dependent rows in (Ag — M). The 
nonsingularity of Ki implies that (ai Ab — M ) has full row rank and there must 
exist a vector v such that v'^ai ^ 0, v'^Ab = 0 and v'^M = 0. If u is scaled so that 
v'^ai = —Azi, then (0,0,—u) must be a solution of (2.17). It follows that Axi = 0, 
V = Ay, and (0, Ay) is an eigenvector of Kb associated with a zero eigenvalue. The 
nonsingularity of Ki implies that v is unique given the value of the scalar Azi, and 
hence the zero eigenvalue has multiplicity one. 

Conversely, Axi = 0 implies that {Axb, Ay) is a null vector for Kb- However, if 


Kb is nonsingular, then the vector is zero, contradicting (2.16a). It follows that Kb 
must be singular. | 


3. A Primal Active-Set Method for Convex QP 

In this section a primal-feasible method for convex QP is formulated. Each iteration 
begins and ends with a point (x, y, z) that satishes the conditions 

Hx + c- A^y - z = 0, Xn + Qn = 0, Xs + gs > 0, 

Ax-\-Aly - b = 0, ZB + rB = 0, 

for appropriate second-order consistent bases. The purpose of the iterations is to 
drive (x, y, z) to optimality by driving the dual variables to feasibility (i.e., by driving 
the negative components of Zn + Vn to zero). Methods for hnding 13 and J\f at the 
initial point are discussed in Section 

An iteration consists of a group of one or more consecutive subiterations during 
which a specihc dual variable is made feasible. The hrst subiteration is called the 
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base subiteration. In some cases only the base subiteration is performed, but, in 
general, additional intermediate subiterations are required. 

At the start of the base subiteration, an index I in the nonbasic set M is identified 
such that zi + ri < 0. The idea is to remove the index I from J\f (i.e., J\f J\f\ {/}) 
and attempt to increase the value of zi + ri by taking a step along a primal-feasible 
direction (Axi, Axb, Ay, Azi). The removal of I from Af implies that BU{1} UAf = 
{1, 2, ..., n} with B second-order consistent. This implies that Kg is nonsingular 


and the (unique) search direction may be computed as in (2.11) with Axi = 1. 

If Azi > 0, the step = —{zi + ri)/Azi gives zi + a^Azi + ri = 0. Otherwise, 
Azi = 0, and there is no finite value of a. that will drive zi P aAzi P ri to its bound, 
and is defined to be -l-oo. Proposition IAT tI of the Appendix implies that the case 


Azi = 0 corresponds to the primal objective function being linear and decreasing 
along the search direction. 

Even if Azi is positive, it is not always possible to take the step and remain 
primal feasible. A positive step in the direction {Axi, Axg, Ay, Azi) must increase 
xi from its bound, but may decrease some of the basic variables. This makes it 
necessary to limit the step to ensure that the primal variables remain feasible. The 
largest step length that maintains primal feasibility is given by 


Q^max — 


mm 

i:Axi<0 


Qi 


-Axi ' 

If «max is finite, this value gives xj- P a^axAxk -|- = 0, where k is the index k = 

argminj.^ 2 ,,<Q {xiP qi) / {—Axi). The overall step length is then a = min (a^,amax)- 
An infinite value of a indicates that the primal problem (PQPg ,,) is unbounded, 
or, equivalently, that the dual problem (DQPg,,) is infeasible. In this case, the 
algorithm is terminated. If the step a = is taken, then zi P aAzi -|- n = 0, the 
subiterations are terminated with no intermediate subiterations and B B L) {1}. 
Otherwise, a = a m ax . and the basic and nonbasic sets are updated as B B\ {k} 
and AA ■<— Mu{k} giving a new partition BU{l}UAf = {1,2,..., n}. In order to show 
that the equations associated with the new partition are well-defined, it is necessary 
to show that allowing z^ to move does not give a singular AT;. Proposition A.5 of 


the Appendix shows that the submatrix Ki associated with the updated B and Af 
is nonsingular for the cases Azi > 0 and Azi = 0. 

Because the removal of k from B does not alter the nonsingularity of Ki, it is 


possible to add I to B and thereby define a unique solution of the system (2.4). 
However, if zi P ri < 0, additional intermediate subiterations are required to drive 
zi P ri to zero. In each of these subiterations, the search direction is computed 
by choosing Azi = 1 in Proposition 2.4 The step length is given by = 


— {ziPri)/Azi as in the base subiteration above, but now a^ is always finite because 
Azi = 1. Similar to the base subiteration, if no constraint is added, then ziPa^^AziP 
r; = 0. Otherwise, the index of another blocking variable k is moved from B to N. 
Proposition |A.5| implies that the updated matrix Ki is nonsingular at the end of an 
intermediate subiteration. As a consequence, the intermediate subiterations may be 
repeated until zi P ri is driven to zero. 

At the end of the base subiteration or after the intermediate subiterations are 
completed, it must hold that zi P ri = 0 and the final Ki is nonsingular. This 
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implies that a new iteration may be initiated with the new basic set BU {/} defining 
a nonsingular K^- 

The primal active-set method is summarized in Algorithm |3.1| below. The con¬ 
vergence properties of Algorithm |3.1| are established in Section which concerns a 


general primal algorithm that includes Algorithm 3.1 as a special case. 


4. A Dual Active-Set Method for Convex QP 

Each iteration of the dual active-set method begins and ends with a point (x, y, z) 
that satisfies the conditions 


Hx + c — A^y — z = 0, 
Ax My — 6 = 0, 


Xn P Qn — 0, 

Zb ~\- f B = 0, 


-Ziv + rjv ^ 0, 


(4.1) 


for appropriate second-order consistent bases. For the dual method, the purpose is 
to drive the primal variables to feasibility (i.e., by driving the negative components 
of X -|- g to zero). 

An iteration begins with a base subiteration in which an index I in the basic 
set B is identified such that xi + qi < 0. The corresponding dual variable zi may 
be increased from its current value zi = —ri by removing the index I from B, and 
defining B B\ {/}. Once I is removed from B, it holds that B U {/} U J\f = {1, 


2,..., n}. The resulting matrix Ki of (2.10) is nonsingular, and the unique direction 


{Axi, Axb, Ay) may be computed with Azi = 1 in Proposition 2.4 


If Axi > 0, the step = —(x; -|- qi)/Axi gives x; -|- a^Axi qi = 0. Otherwise, 


Axi = 0 and Proposition |A.7| of the Appendix implies that the dual objective 
function is linear and increasing along {Ax,Ay,Az). In this case = -|-cx). As 
xi -\- qi is increased towards zero, some nonbasic dual variables may decrease and the 
step must be limited by Umax = (zi -f ri)(—Azi) to maintain feasibility 

of the nonbasic dual variables. This gives the step a = min (a^, amax)- If a = -|-oo, 
the dual problem is unbounded and the iteration is terminated. This is equivalent 
to the primal problem (PQPg ,,) being infeasible. If a = a^, then xi + aAxi + qi = 0. 
Otherwise, it must hold that a = Umax and J\f and B are redefined as J\f = J\f \ {A:} 
and B = Be {k}, where k is the index k = argmin^.^^^^Q (zi + ri)j{—Azi). The 
partition at the new point satisfies B U {/} U A/" = {1, 2, ..., n}. Proposition A.6 of 
the Appendix shows that the new Kb is nonsingular for both of the cases Axi > 0 
and Ax I = 0. 

If XI qi < 0 at the new point, then at least one intermediate subiteration 
is necessary to drive x; -|- qi to zero. The nonsingularity of Kb implies that the 
search direction may be computed with Axi = 1 in Proposition |2.3[ As in the 
base subiteration, the step length is = —{xi + qi)/Axi, but in this case can 
never be infinite because Axi = 1. If no constraint index is added to B, then 
xi -|- aAxi qi = 0. Otherwise, the index A; of a blocking variable is moved from M 
to B. Proposition |A.6| of the Appendix implies that the updated Kb is nonsingular 
at the end of an intermediate subiteration. Once x; -|- qi is driven to zero, the index 
I is moved to Af and a new iteration is started. 
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Algorithm 3.1 A primal active-set method for convex QP. 


Find {x,y,z) satisfying conditions (3.1) for some second-order consistent basis B] 
^vhile 3 I : zi + ri < 0 do 
Af Af\ {Q; 

primal_base(;B, M, I, X, y, z); 
while zi + ri < 0 do 

primal_intermediate(.B, Af, I, X, y, z); 

end while 
B^BC{1}- 
end while 


[returns B, Af, x, y, z] 
[returns B, Af, x, y, z] 


function primal_base(.B, Af, I, x, y. 

Solve 


v) 

-^y J 


Azm h^^Axl + Hl^Axg - A^Ay, 
Azi <— hiiAxi + hffgiAxg — af Ay; 


a* < - {zi + ri)/Azi; 

Omax ^ min {Xi + qi)/{-Axi); /c ^ argmin 
i-.Axi<0 i:Axi<0 


hBi \ _ 


[Azi > 0] 
[«* —1-00 if Azi = 0] 
[xi qi)/{-Axi); 


a min (a*,amax); 
if a = -|-oo then 

stop; [(DQPq,,) is infeasible] 

end if 

xi xi + aAxi; Xg Xg + aAxg] 

y^y + aAy; zi ^ zi + aAzi; z^ ^ + aAzr^; 

if z; -|- n < 0 then 


B^B\{k}; Af^Nc{k}; 


end if 

return B, Af, x, y, z; 

end function 


function PRIMAL_INTERMEDIATE(i3, Af, I, X, y, z) 

(K af \ / Axi\ / 1 \ 

Azi ^ 1; Solve Hgg Af Axg = 0 ; [Axi > 0] 

\ai As -mJ \-Ay ) \oJ 

Azn ^ HnA^i + Hf^Axg - AfAy; 
a* ^- {zi + ri); 

Omax ^ , min {xi qf)/{-Axi); k ^ argmin (x* qi)/{-Axi); 
i:Axi<0 i:Axi<0 

ex i min , Omax) : 

xi xi + aAxi; Xg Xg + aAxg] 

y^y + aAy; zi^zi + aAzi; ZM^z^ + aAzu; 

if z; -|- r; < 0 then 

B^B\{k]; N^Nc{k]; 
end if 

return B, Af, x, y, z; 

end function 
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The dual active-set method is summarized in Algorithm 4.1 below. Its conver¬ 


gence properties are discussed in Section 5.5 


5. Combining Primal and Dual Active-Set Methods 


The primal active-set method proposed in Sectionmay be used to solve (PQPg,,) 
for a given initial second-order consistent basis satisfying the conditions g. An 
appropriate initial point may be found by solving a conventional phase-1 linear 
program. Alternatively, the dual active-set method of Section may be used in 
conjunction with an appropriate phase-1 procedure to solve the quadratic program 
(PQPq ,.) for a given initial second-order consistent basis satisfying the conditions 
(4.1). In this section a method is proposed that provides an alternative to the 


conventional phase-1/phase-2 approach. It is shown that a pair of coupled quadratic 
programs may be created from the original by simultaneously shifting the bound 
constraints. Any second-order consistent basis can be made optimal for such a 
primal-dual pair of shifted problems. The shifts are then updated using the solution 
of either the primal or the dual shifted problem. An obvious application of this 
approach is to solve a shifted dual QP to dehne an initial feasible point for the 
primal, or vice-versa. This strategy provides an alternative to the conventional 
phase-1/phase-2 approach that utilizes the QP objective function while Ending a 
feasible point. 


5.1. Finding an initial second-order-consistent basis 


For the methods described in Section 5.2 below, it is possible to dehne a simple 
procedure for finding the initial second-order consistent basis B such that the ma¬ 
trix Kg of (2.6) is nonsingular. The required basis may be obtained by Ending a 
symmetric permutation 77 of the “full” KKT matrix K such that 


77^7777 = 77 '^ 


77 

A -M 


77 = 


Hgg 

Al 

^BN 

Ag 

-M 

Ajv 

ttT 

BN 

■^N 

^NN 


(5.1) 



Such methods use symmetric interchanges that implicitly form the nonsingular ma¬ 
trix Kg by deferring singular pivots. In this case, Kg may be deEned as any sub¬ 
matrix of the largest nonsingular principal submatrix obtained by the factorization. 
(There may be further permutations within 77 that are not relevant to this dis- 
for further details, see, e.g.. 


20,21 28,^.) The permutation 77 deEnes 


cussion; 


the initial B-N partition of the columns of A, i.e., it deEnes an initial second-order 
consistent basis. 
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Algorithm 4.1 A dual active-set method for convex QP. 


Find {x,y,z) satisfying conditions (4.1) for some second-order consistent basis B] 
while 3 I : xi + qi < 0 do 
B^B\{iy, 

dual_base(;B, a/", I, X, y, z); [Base subiteration] 

while xi + qi < 0 do 


dual_intermediate(.B, J\f, I, X, y, z); 

end while 

J\f^Afc{l}-, 

end while 


[Intermediate subiteration] 


function dual_base(.B, Af, I, x, y, z) 

hu hi af\/ Axi\ /1\ 

Azi ^ 1; Solve hgi Hbb = 0 h - 0] 

\ ai Ab -m) \ -Ay ) \0/ 

Az]\i i h^^Axj^ PIAxB 

a* < - {xi + qi)/Axi] [a^ -Foo if Axi = 0] 

Omax ^ min {Zi + Vi)/{-Azi)\ k ^ argmin {zi + ri)/{-Azi)] 

r.Azi <0 r.Azi <0 

a min (a*,amax); 

if a = -t-oo then 

stop; [(PQPq,,) is infeasible] 

end if 

xi xi + aAxi; Xb Xb + uAxb] 

y-^y + aAy; zizi + aAzf, z^z^ + aAz^', 

if xi + qi < 0 then 

B^BU{ky, M^M\{ky 

end if 

return B, Af, x, y, z] 

end function 


function DUAL_INTERMEDIATE(i3, A/, I, X, y, z) 


Ax I ■<— 1; Solve 


Hbb 

Ab 


Al 

-M 


Axb 

-Ay 


Azi ■<— hiiAxi + HIAxb — af Ay; 

Azjs! i hf^^Axi P HBfjAxB Aj^Ay, 

a* ^- {xi + qi); 

Omax ^ min {zi + Vi)/ {-Azi); k ^ argmin 

r.Azi <0 


a 


mm I 


hBi \ ^ 
ai)' 


(zi + ri)/{-Azi); 


xi xi + aAxi; Xb Xb + oAxb] 
y^y + aAy; zi^zi + aAzi; ZM^z^ + aAzu; 
xi + qi < then 

B^BC{ky Af^Af\{ky 


end if 

return B, Af, x, y, z; 

end function 


[Azi > 0] 
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5.2. Initializing the shifts 

Given a second-order consistent basis, it is straightforward to create shifts 
and corresponding {x, y, z) so that > 0, > 0 and (x, y, z) are optimal for 

(PQPg(o) ,,( 0 )) and (DQPg(o) ,,(o)). First, choose nonnegative vectors and 
(Obvious choices are = 0 and = 0.) Define Zg = —r^g\ and 


(2.7). Finally, let > max{—Xs,0} and > max{—Zjv,0}. Then, it follows 
that X, y and z are optimal for the problems (PQP (o) ^(o)) and 


solve the nonsingular KKT-system (2.5) to obtain Xg and y, and compute Zpj from 


(0) 


2.1 


from Proposition 

(DQPg(o) r(o)), with q^^^ > 0 and > 0. If q^^^ and are zero, then x, y and z 
are optimal for the original problem. 


5.3. Solving the original problem by removing the shifts 

The original problem may now be solved as a pair of shifted quadratic programs. 
Two alternative strategies are proposed. The first is a “primal first” strategy in 
which a shifted primal quadratic program is solved, followed by a dual. The second 
is an analogous “dual first” strategy. 

The “primal-first” strategy is summarized as follows. 

(0) Find B, J\f, q^^\ x, y, z, as described in Sections 


5.1 


and 


5.2 


(1) Set = 0. Solve (PQPq,o) using the primal active-set method. 

(2) Set q^^'^ = 0, = 0. Solve (DQPq^o) using the dual active-set method. 

In steps (1) and (2), the initial B-Af partition and initial values of x, y, and z are 
defined as the final B-M partition and final values of x, y, and z from the preceding 
step. 

The “dual-first” strategy is defined in an analogous way. 

(0) Find B, J\f, q^^\ x, y, z, as described in Section 


5.1 


and 


5.2 


(1) Set = 0, Solve (DQPg,,) using the dual active-set method. 

(2) Set q^‘^'i = 0, = 0. Solve (PQPq^o) using the primal active-set method. 

As in the “primal-first” strategy, the initial B-M partition and initial values of x, y, 
and z for steps (1) and (2), are defined as the final B-M partition and final values 
of X, y, and 2; from the preceding step. 

(The strategies of solving two consecutive quadratic programs may be generalized 
to a sequence of more than two quadratic programs, where we alternate between 
primal and dual active-set methods, and eliminate the shifts in more than two steps.) 
In order for these approaches to be well-defined, a simple generalization of the 


primal and dual active-set methods of Algorithms 3.1 and |4.1|is required. 


5.4. Relaxed initial conditions for the primal QP method. 

For A lgorithm [3T| the initial values of B, M, q, r, x, y, and 2 ; must satisfy conditions 
(3.1). However, the choice of r = = 0 in Step (2) of the dual-first strategy 

may give some negative components in the vector Zg rg. This possibility may 
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be handled by defining a simple generalization of Algorithm 3.1 that allows initial 
points satisfying the conditions 


Hx + c — A^y — z = 0, 
Ax + My — b = 0, 


xjv + giv = 0, 

Zb + rs < 0, 


Xb PQb > 0, 


(5.2) 


instead of the conditions (3.1). In Algorithm 3.1, the index I identified at the start 


of the primal base subiteration is selected from the set of nonbasic indices such that 
Zj + rj < 0. In the generalized algorithm, the set of eligible indices for I is extended 
to include indices associated with negative values of Zb+Vb- If the index I is deleted 
from B, the associated matrix Ki is nonsingular, and intermediate subiterations are 
executed until the updated value satisfies zi + ri = 0. At this point, the index I is 


returned B. The method is summarized in Algorithm 5.1 


Algorithm 5.1 A primal active-set method for convex QP. 


Find (x,y,z) satisfying conditions (5.2) for some second-order consistent basis B] 
while 3 I : zi + ri < 0 do 

if I £ J\f then 

Af Af\ {1}', 


primal_base(;B, Af, I, X, y, z); 

else 

B^B\{1}; 

end if 

while zi + ri < 0 do 

primal_intermediate(.B, Af, I, X, y, z) 

end while 
B^BU{1}; 
end while 


[returns B, Af, x, y, zj 


[returns B, Af, x, y, zj 


This section concludes with a convergence result for the primal method of Al¬ 
gorithm 5.1 In particular, it is shown that the algorithm is well-defined, and ter¬ 
minates in a finite number of iterations if (PQP^ ,,) is nondegenerate. We define 
nondegeneracy to mean that a nonzero step in the x-variables is taken at each iter¬ 


ation of Algorithm 5.1 that involves a base subiteration. A sufficient condition on 
(PQPq,.) for this to hold is that the gradients of the equality constraints and active 
bound constraints are linearly independent at each iterate. See, e.g., Fletcher pQ 


for further discussion of these issues. As the active-set strategy uses the same crite¬ 
ria for adding and deleting variables as those used in the simplex method, standard 
pivot selection rules used to avoid cycling in linear programming, such as lexico¬ 
graphical ordering, least-index selection or perturbation may be applied directly to 
the method proposed here (see, e.g.. Bland [^, Charnes [^, Dantzig, Orden and 
Wolfe flTl, and Harris |4^). 


Theorem 5.1. Given a primal-feasible point (xo,yo,zo) satisfying conditions (5.2) 
for a second-order consistent basis Bq, then Algorithm |5.1| generates a sequence of 
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second-order consistent bases {Bj}j^o. Moreover, if problem (PQPg,,) is nondegen¬ 
erate, then Algorithm 5.1 finds a solution o/(PQP^,.) or determines that (DQPg ,,) 


is infeasible in a finite number of iterations. 


Proof. Assume that {x,y,z) satisfies the conditions (5.2) for the second-order 
consistent basis B. Propositions |2.3| and |2.4| imply that the KKT matrices associated 
with subsequent base and intermediate iterations are nonsingular, in which case each 
basis is second-order consistent. Let B^ denote the index set B^ = {i ^ B : Zi-\- Ti <C 
0}, and let r be the vector Fj = ri, i ^ B^, and Fj = —Zi, i G B^. These definitions 
imply that Fj = —Zi > —Zi Zi ri = ri, for every i £ B^. It follows that r > r, 
and the feasible region of (DQP ) is a subset of the feasible region of (DQP jr). In 


addition, if r is replaced by r in (3.1), the only difference is that Zg + rg = 0, i.e.. 


the initial point for (5.2) is a stationary point with respect to (PQP 


q,r/ 


The first step of the proof is to show that after a finite number of iterations of 


Algorithm 5.1, one of three possible events must occur: (i) the cardinality of the 
set B^ is decreased by at least one; (ii) a solution of problem (PQP^,,) is found; or 
(hi) (DQPq,,) is declared infeasible. The proof will also establish that if (i) does not 
occur, then either (ii) or (hi) must hold after a finite number of iterations. 

Assume that (i) never occurs. This implies that the index I selected in the base 
subiteration can never be an index in B^ because at the end of such an iteration, it 
would belong to B with zi-\-ri = 0, contradicting the assumption that the cardinality 
of B"^ never decreases. For the same reason, it must hold that k ^ B^ for every index 
k selected to be moved from .6 to AA in any subiteration, because an index can only 
be moved from AA to .6 by being selected in the base subiteration. These arguments 
imply that Zi = —ri, with i £ B^, throughout the iterations. It follows that the 
iterates may be interpreted as being members of a sequence constructed for solving 
(PQPg jr) with a fixed F, where the initial stationary point is given, and each iteration 
gives a new stationary point. The nondegeneracy assumption implies that a Ax 0 
for at least one subiteration. For the base subiteration, Axi > 0, and it follows from 
Proposition |2.4| that Ax 7 ^ 0 if and only if Axi > 0 for an intermediate subiteration. 
Therefore, Proposition A.7 shows that the objective value of (PQPg jr) is strictly 


decreasing for a subiteration where aAx 0. In addition, the objective value of 
(PQPg jr) is nonincreasing at each subiteration, so a strict overall improvement of the 
objective value of (PQPg^r) is obtained at each iteration. As there are only a finite 

either solves (PQP^ 7 


5.1 


number of stationary points. Algorithm 

(DQPq jr) is infeasible after a finite number of iterations. If (PQP 


or concludes that 
is solved, then 


Zn P "<'n > 0, because Fj = for f £ M. Hence, Algorithm |5.1| can not proceed 
further by selecting an I £ J\f, and the only way to reduce the objective is to select 
an Z in H such that Zj -£ rj < 0. Under the assumption that (i) does not occur, it 


= 0. However, in this case (PQP^,.) 

declares (DQP^ jr) 


5.1 


must hold that no eligible indices exist and B^ 
has been solved with r = r, and (ii) must hold. If Algorithm 
to be infeasible, then (DQP^,,) must also be infeasible because the feasible region 
of (DQPg,.) is contained in the feasible region of (DQP^ jr). In this case (DQP^,,) is 
infeasible and (hi) occurs. 

Finally, if (i) occurs, there is an iteration at which the cardinality of B^ decreases 
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and an index is removed from . There may be more than one such index, but 
there is at least one I moved from B^ to B\B^, or one k moved from B^ to Af. 
In either case, the cardinality of B^ is decreased by at least one. After such an 
iteration, the argument given above may be repeated for the new set B^ and new 
shift r. Applying this argument repeatedly gives the result that the situation (i) 
can occur only a hnite number of times. 

It follows that (ii) or (iii) must occur after a finite number of iterations, which 
is the required result. | 


5.5. Relaxed initial conditions for the dnal QP method. 

Analogous to the primal case, the choice of q = = 0 in Step (2) of the primal- 

first strategy may give some negative components in the vector Xjv + ^at. In this 


case, the conditions (4.1) on the initial values of B, Af, q, r, x, y, and z are relaxed 
so that 

Hx + c — A^y — z = 0, Xpf + qj^ < 0, 

Ax -|- My — 6 = 0, 2 b + Ts = 0, Pr^ >2^. 


(5.3) 


Similarly, the set of eligible indices may be extended to include indices associated 
with negative values of x^ + qN- If the index I is from Af, the associated matrix 
Kg is nonsingular, and intermediate subiterations are executed until the updated 
value satishes xi + qi = 0. At this point, the index I is returned Af. The method is 


summarized in Algorithm 5.2 


Algorithm 5.2 A dual active-set method for convex QP. 


Find {x,y,z) satisfying conditions (5.3) for some second-order consistent B] 
while 3 I : xi P qi < 0 do 


if I € B then 

B^B\{iy, 

dual_base(;B, Af, l, X, y, z); [Base subiteration] 

else 


Af Af\ {Z}; 

end if 

while xi P qi < 0 do 

dual_intermediate(.B, Af, I, X, y, z); [Intermediate subiteration] 

end while 
Af^Afc{l}-, 
end while 


A convergence result analogous to Theorem 5.1 holds for the dual algorithm. In 
this case, the nondegeneracy assumption concerns the linear independence of the 
gradients of the equality constraints and active bounds for (DQPg ,,). 


Theorem 5.2. Given a dual-feasible point {xo,yo,zo) satisfying conditions (5.3) 
for a second-order consistent basis Bq, then Algorithm 5^ generates a sequence 
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of second-order consistent bases {i3j}j>o. 
degenerate, then Algorithm 


Moreover, if problem (DQPg^) is non 


5.2 


either solves (DQP„ or concludes that (PQP„ r) 


infeasible in a finite number of iterations. 


Proof. The proof mirrors that of Theorem 5.1 for the primal method. 


6. Practical Issues 

As stated, the primal quadratic program has lower bound zero on the x-variables. 
This is for notational convenience. This form may be generalized in a straightfor¬ 
ward manner to a form where the x-variables has both lower and upper bounds on 
the primal variables, i.e., < x < bu, where components of 6^ can be —oo and 

components of bu can be + 00 . Given primal shifts qu and Qu, and dual shifts and 
ru, we have the primal-dual pair 

minimize Ix^Hx \iP"My -\- c^x (ru — ru)'^x 

(PQPg^^) 

subject to Ax + My = 6, bu — Qu < x < bu + Qu, 


and 


maximize -\x^Hx - ly^My -[- b^y -\- (bu - quVzu - (bu + quVzu 

(DQPq_^) 

subject to —Hx + A^y Zu — Zu = c, Zu > —ru, Zu > —ru- 


An infinite bound has neither a shift nor a corresponding dual variable. For example, 
if the jth components of bu and bu are infinite, then the corresponding variable Xj 
is free. In the procedure given in Section 5.1 for finding the first second-order 
consistent basis B, it is assumed that variables with indices not selected for B are 
initialized at one of their bounds. As a free variable has no finite bounds, any index 
j associated with a free variable should be selected for B. However, this cannot be 
guaranteed in practice, and in the next section it is shown that the primal and dual 
QP methods may be extended to allow a free variable to be fixed temporarily at 
some value. 

If the QP is defined in the general problem format of Section then any free 
variable not selected for B has no upper or lower bound and must be temporarily 
fixed at some value xj = xj (say). The treatment of such “temporary bounds” in¬ 


volves some additional modifications to the primal and dual methods of Sections 5.4 
and 15.51 


Each temporary bound Xj = Xj defines an associated dual variable Zj with initial 
value Zj. As the bound is temporary, it is treated as an equality constraint, and 
the desired value of Zj is zero. Initially, an index j corresponding to a temporary 
bound is assigned a primal shift qj = 0 and a dual shift rj = —zj, making Xj and zj 
feasible for the shifted problem. In both the primal-first and dual-first approaches, 
the idea is to drive the z^-variables associated with temporary bounds to zero in the 
primal and leave them unchanged in the dual. 
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In a primal problem, regardless of whether it is solved before or after the dual 


problem, an index j corresponding to a temporary bound for which Zj 7 ^ 0 is consid¬ 
ered eligible for selection as I in the base subiteration, i.e., the index can be selected 


regardless of the sign of Zj. Once selected, Zj is driven to zero and j belongs to B 
after such an iteration. In addition, as xj has no finite bounds, j will remain in B 
throughout the iterations. Hence, at termination of a primal problem, any index j 
corresponding to a temporarily bounded variable must have Zj = 0. If the maximum 
step length at a base subiteration is infinite, the dual problem is infeasible, as in the 
case of a regular bound. 

In a dual problem, the dual method is modified so that the dual variables as¬ 
sociated with temporary bounds remain fixed throughout the iterations. At any 
subiteration, if it holds that Azj / 0 for some temporary bound, then no step is 
taken and one such index j is moved from M to B. Consequently, a move is made 
only if Azj = 0 for every temporary bound j. It follows that the dual variables for 
the temporary bounds will remain unaltered throughout the dual iterations. Note 
that an index j corresponding to a temporary bound is moved from AA to H at most 
once, and is never moved back because the corresponding Xj-variable has no finite 
bounds. If the maximum step length at a base subiteration is infinite, it must hold 
that Azj = 0 for all temporary bounds j, and the primal problem is infeasible. 

The discussion above implies that a pair of primal and dual problems solved 
consecutively will terminate with Zj = 0 for all indices j associated with temporary 
bounds. This is because Zj is unchanged in the dual problem and driven to zero in 
the primal problem. 

7. Numerical Examples 

This section concerns a particular formulation of the combined primal-dual method 
of Section in which either a “primal-hrst” or “dual-first” strategy is selected based 
on the initial point. In particular, if the point is dual feasible, then the “dual-first” 
strategy is used, otherwise, the “primal-first” strategy is selected. Some numerical 
experiments are presented for a simple Matlab implementation applied to a set 
of convex problems from the CUTEst test collection (see Bongartz, et al. |^, and 
Gould, Orban and Toint (4^[4^ ). 

7.1. The test problems 

Each QP problem in the CUTEst test set may be written in the form 



minimize Ix'^Hx -|- c'^x subject to 


where I and u are constant vectors of lower and upper bounds, and A has dimension 
mx n. In this format, a fixed variable or equality constraint has the same value for 
its upper and lower bound. Each problem was converted to the equivalent form 


minimize \x'^Hx + (Fx subject to Ax — s = 0, 

rr> a ^ 



(7.1) 
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where s is a vector of slack variables. With this formulation, the QP problem involves 
simple upper and lower bounds instead of nonnegativity constraints. It follows that 
the matrix M is zero, but the full row-rank assumption on the constraint matrix is 
satisfied because the constraint matrix A takes the form (^A — l) and has rank m. 

Numerical results were obtained for a set of 121 convex QPs in standard interface 
format (SIF). The problems were selected based on the dimension of the constraint 


matrix A in (7.1). In particular, the test set includes all QP problems for which the 


smaller of m and n is of the order of 500 or less. This gave 121 QPs ranging in size 
from BQPIVAR (one variable and one constraint) to LINCONT (1257 variables and 419 
constraints). 


7.2. The implementation 

The combined primal-dual active-set method was implemented in Matlab as Al¬ 
gorithm PDQP. For illustrative purposes, results were obtained for PDQP and the QP 
solver SQOPT 35 , which is a Fortran implementation of a conventional two-phase 


(primal) active-set method for large-scale QP. Both PDQP and SQOPT use the method 


of variable reduction, which implicitly transforms a KKT system of the form (2.5) 


into a block-triangular system. The general QP constraints Ax — s = 0 are parti¬ 
tioned into the form Bx^ -|- Sx^ AffX^f = 0, where B is square and nonsingular, 
with Ag = {B S) and Xg = (x^,x®). The vectors x®, x^ are the associated 
basic, superbasic, and nonbasic components of (x, s) (see Gill, Murray and Saun¬ 
ders 34 ). If 77 denotes the Hessian H of (7.1) augmented by the zero rows and 


columns corresponding to the slack variables, then the reduced Hessian Z^HZ is 
defined in terms of the matrix Z such that 



where P permutes the columns of (A ~ ™to the order (H S Aj^ ). The 
matrix Z is used only as an operator, i.e., it is not stored explicitly. Products of the 
form Zv or Z'^u are obtained by solving with B or B^. With these definitions, the 
resulting block lower-triangular system has diagonal blocks Z^HZ, B and B^. 

The initial nonsingular B is identified using an LU factorization of . The 
resulting Z is used to form Z^HZ, and a partial Cholesky factorization with inter¬ 
changes is be used to find an upper-triangular matrix R that is the factor of the 
largest nonsingular leading submatrix of Z^HZ. If Zg denotes the columns of Z 
corresponding to R, and Z is partitioned as Z = i^Zg Z^'j, then the index set B 
consisting of the union of the column indices of B and the indices corresponding to 
Zg defines an appropriate initial second-order consistent basis. 

All SQOPT runs were made using the default parameter options. Both PDQP and 
SQOPT are terminated at a point (x, y, z) that satishes the optimality conditions of 
Proposition |2 .1 1 modified to conform to the constraint format of (7.1). The feasibility 


and optimality tolerances are given by Cfea = 10 ° and Copt = 10 °, respectively. 
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For a given eopt, PDQP and SQOPT terminate when 

> -eoptllylloo if Xi > -ii, i G Af] 

Zi < Coptllylloo if Xi < Ui,ie M. 

Both PDQP and SQOPT use the EXPAND anti-cycling procedure of Gill et al. to 
allow the variables (x, s) to move outside their bounds by as much as Cfea- The 
EXPAND procedure does not guarantee that cycling will never occur (see Hall and 
McKinnon for an example). Nevertheless, in many years of use, the authors 
have never known EXPAND to cycle on a practical problem. 


7.3. Numerical results 


PDQP and SQOPT were applied to the 121 problems considered in Section 7.1 A 


summary of the results is given in Table The first four columns give the name of 
the problem, the number of linear constraints m, the number of variables n, and the 
optimal objective value Objective. The next two columns summarize the SQOPT 
result for the given problem, with Phsl and Itn giving the phase-one iterations 
and iteration total, respectively. The last four columns summarize the results for 
PDQP. The first column gives the total number of primal and dual iterations Itn. The 
second column gives the order in which the primal and dual algorithms were applied, 
with PD indicating the “primal-first” strategy, and DP the “dual-first” strategy. The 
hnal two columns, headed by p-Itn, and d-Itn, give the iterations required for the 
primal method and the dual method, respectively. 

Of the 121 problems tested, two (LINCDNT and NASH) are known to be infeasible. 
This infeasibility was identified correctly by both SQOPT and PDQP. In total, SQOPT 
solved 117 of the remaining 119 problems, but declared (incorrectly) that problems 
RDW2D51U and RDW2D52U are unbounded. PDQP solved the same number of problems, 
but failed to achieve the required accuracy for the problems RDW2D51B and RDW2D52F. 
In these two cases, the final objective values computed by PDQP were 1.0947648E-02 
and 1.0491239E-02 respectively, instead of the optimal values 1.0947332e-02 and 
1.0490828e-02. (The five RDW2D5* problems in the test set are known to be difficult 
to solve, see Gill and Wong [40].) 

Figure gives a performance profile (in log2 scale) for the iterations required by 
PDQP and SQOPT. (For more details on the use of performance profiles, see Dolan and 
More |19|.) The hgure profiles the total iterations for PDQP, the number of phase- 
2 iterations for SQOPT, and the sum of phase-1 and phase-2 iterations for SQOPT. 
Some care must be taken when interpreting the results in the prohle. First, the 
CUTEst test set contains several groups made up of similar variants of the same 
problem. In this situation, the profiles can be skewed by the fact that a method will 
tend to exhibit similar performance on all the problems in the group. For example, 
PDQP performs significantly better than SQOPT on all four JMLBRNG* problems, but 
signihcantly worse on all 12 LISWET* problems. 

Second, the phase-1 search direction for SQOPT requires the computation of the 
vector —ZZ'^g{x), where g{x) is the gradient of the sum of infeasibilities of the 
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bound constraints at x. This implies that a phase-1 iteration for SQOPT requires 
solves with B and , compared to solves with B, B"^ and Z^HZ for a phase-2 
iteration. As every iteration for PDQP requires the solution of a KKT system, if the 
number of superbasic variables is not small, a phase-1 iteration of SQOPT requires 
considerably less work than an iteration of PDQP. It follows that the total iterations 
for PDQP and SQOPT are not entirely comparable. In particular a profile that would 
provide an accurate comparison with PDQP lies somewhere in-between the two SQOPT 
profiles shown. 

Notwithstanding these remarks, the profile indicates that PDQP has comparable 
overall performance to a primal method that ignores the objective function while 
finding an initial feasible point. This provides some preliminary evidence that a 
combined primal-dual active set method can be an efficient and reliable alternative 
to conventional two-phase active-set methods. The relative performance of the pro¬ 
posed method is likely to increase when solving a sequence of related QPs for which 
the initial point for one QP is close to being the solution for the next. In this case, 
regardless of whether a primal or dual method is being used to solve the QP, the 
initial point may start off being primal or dual feasible, or the number of primal 
or dual infeasibilities may be small. This is typically the case for QP subproblems 
arising in sequential quadratic programming methods or mixed-integer QP. 

Figure provides a bar graph of the so-called “outperforming factors” for it¬ 
erations, as proposed by Morales 53 . On the x-axis, each bar corresponds to a 


particular test problem, with the problems listed in the order of Table The y-axis 
indicates the factor (log 2 scaled) by which one solver outperformed the other. A 
bar in the positive region indicates that PDQP outperformed SQOPT. A negative bar 
means SQOPT performed better. A positive (negative) dark grey bar denotes a failure 
in SQOPT (PDQP). Light grey bars denote a zero iteration count for a solver. 


Table 1: Results for PDQP and SQOPT on 121 CUTEst QPs. 


Name 

m 

n 

Objective 

SQOPT 

Phsl Itn 

Itn 

PDQP 

Order P-Itn 

D-Itn 

ALLINQP 

50 

100 

-9.1592833E+00 

0 

45 

65 

PD 

63 

2 

AUG2DqP 

100 

220 

1.7797215E+02 

8 

116 

440 

PD 

326 

114 

AUG3D 

27 

156 

8.3333333E-02 

0 

45 

45 

DP 

0 

45 

AVGASA 

10 

8 

-4.6319255E+00 

5 

8 

5 

DP 

0 

5 

AVGASB 

10 

8 

-4.4832193E+00 

5 

8 

7 

DP 

0 

7 

BIGGSBl 

1 

100 

1.5000000E-02 

0 

103 

101 

PD 

101 

0 

BQPIVAR 

1 

1 

O.OOOOOOOE+00 

0 

1 

1 

DP 

0 

1 

BQPGABIM 

1 

50 

-3.7903432E-05 

0 

36 

7 

PD 

7 

0 

BQPGASIM 

1 

50 

-5.5198140E-05 

0 

40 

8 

PD 

8 

0 

CHENHARK 

1 

100 

-2.0000000E+00 

0 

132 

32 

DP 

0 

32 

CVXBQPl 

1 

100 

2.2725000E+02 

0 

100 

119 

DP 

2 

117 

cvxqpi 

50 

100 

1.1590718E+04 

5 

67 

91 

DP 

1 

90 

cvxqp2 

25 

100 

8.1209404E+03 

2 

82 

85 

DP 

2 

83 

cvxqp3 

75 

100 

1.1943432E+04 

17 

46 

113 

DP 

2 

111 

DEGENqP 

1005 

10 

O.OOOOOOOE+00 

0 

6 

18 

PD 

18 

0 

DT0C3 

18 

29 

2.2459038E+02 

1 

10 

17 

DP 

0 

17 

DUALl 

1 

85 

3.5012967E-02 

0 

88 

88 

PD 

88 

0 

DUAL2 

1 

96 

3.3733671E-02 

0 

99 

99 

PD 

99 

0 

DUAL3 

1 

111 

1.3575583E-01 

0 

106 

106 

PD 

106 

0 

DUAL4 

1 

75 

7.4609064E-01 

0 

61 

61 

PD 

61 

0 

DUALCl 

215 

9 

6.1552516E+03 

1 

9 

4 

DP 

0 

4 
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Table 1: Results for PDQP and SQOPT on 121 CUTEst QPs. (continued) 


Name 

m 

n 

Objective 

SQOPT 

Phsl Itn 

Itn 

PDQP 

Order P-Itn 

D-Itn 

DUALC2 

229 

7 

3.5513063E+03 

2 

4 

4 

DP 

0 

4 

DUALC5 

278 

8 

4.2723256E+02 

1 

7 

6 

DP 

0 

6 

DUALC8 

503 

8 

1.8309361E+04 

4 

6 

8 

DP 

0 

8 

GENHS28 

8 

10 

9.2717369E-01 

0 

3 

5 

DP 

0 

5 

GMNCASE2 

1050 

175 

-9.9444495E-01 

18 

99 

91 

DP 

0 

91 

GMNCASE3 

1050 

175 

1.5251466E+00 

31 

100 

86 

DP 

0 

86 

GMNCASE4 

350 

175 

5.9468849E+03 

74 

171 

175 

DP 

0 

175 

G0ULDqP2 

199 

399 

9.0045697E-06 

0 

213 

419 

DP 

0 

419 

GOULDQPS 

199 

399 

5.6732908E-02 

0 

200 

406 

PD 

205 

201 

GRIDNETA 

100 

180 

9.5242163E+01 

5 

35 

134 

PD 

81 

53 

GRIDNETB 

100 

180 

4.7268237E+01 

0 

81 

97 

DP 

0 

97 

GRIDNETC 

100 

180 

4.8352347E+01 

6 

93 

153 

DP 

0 

153 

HS3 

1 

2 

O.OOOOOOOE+00 

0 

2 

1 

PD 

1 

0 

HS3M0D 

1 

2 

1.2325951E-32 

0 

2 

1 

PD 

1 

0 

HS21 

1 

2 

-9.9960000E+01 

0 

1 

0 

PD 

0 

0 

HS28 

1 

3 

1.2325951E-32 

0 

2 

0 

PD 

0 

0 

HS35 

1 

3 

l.lllllllE-01 

0 

5 

1 

DP 

0 

1 

HS35I 

1 

3 

l.lllllllE-01 

0 

5 

1 

DP 

0 

1 

HS35M0D 

1 

3 

2.5000000E-01 

0 

1 

0 

PD 

0 

0 

HS44 

6 

4 

-1.5000000E+01 

0 

2 

4 

PD 

4 

0 

HS44NEW 

6 

4 

-1.5000000E+01 

0 

4 

9 

PD 

9 

0 

HS51 

3 

5 

-8.8817841E-16 

0 

2 

0 

DP 

0 

0 

HS52 

3 

5 

5.3266475E+00 

0 

2 

1 

DP 

0 

1 

HS53 

3 

5 

4.0930232E+00 

0 

2 

1 

DP 

0 

1 

HS76 

3 

4 

-4.6818181E+00 

0 

4 

4 

DP 

0 

4 

HS76I 

3 

4 

-4.6818181E+00 

0 

4 

4 

DP 

0 

4 

HS118 

17 

15 

6.6482045E+02 

0 

21 

23 

DP 

0 

23 

HS268 

5 

5 

7.2759576E-12 

0 

8 

0 

PD 

0 

0 

HUES-MQD 

2 

100 

3.4829823E+07 

1 

103 

7 

DP 

0 

7 

HUES!IS 

2 

100 

3.4829823E+09 

1 

103 

7 

DP 

0 

7 

JNLBRNGl 

1 

529 

-1.8004556E-01 

0 

292 

82 

PD 

82 

0 

JNLBRNG2 

1 

529 

-4.1023852E+00 

0 

252 

42 

PD 

42 

0 

JNLBRNGA 

1 

529 

-3.0795806E-01 

0 

292 

292 

PD 

292 

0 

JNLBRNGB 

1 

529 

-6.5067871E+00 

0 

247 

247 

PD 

247 

0 

KSIP 

1001 

20 

5.7579792E-01 

0 

2847 

36 

DP 

0 

36 

LINCONT 

419 

1257 

infeasible 

138 

138* 

304* 

DP 

0 

304 

LISWETl 

100 

106 

2.6072632E-01 

0 

52 

401 

DP 

0 

401 

LISWET2 

100 

106 

2.5876398E-01 

0 

63 

378 

DP 

0 

378 

LISWET3 

100 

106 

2.5876398E-01 

0 

64 

378 

DP 

0 

378 

LISWET4 

100 

106 

2.5876399E-01 

0 

61 

378 

DP 

0 

378 

LISWET5 

100 

106 

2.5876410E-01 

0 

58 

378 

DP 

0 

378 

LISWET6 

100 

106 

2.5876390E-01 

0 

67 

378 

DP 

0 

378 

LISWET7 

100 

106 

2.5895785E-01 

0 

68 

378 

DP 

0 

378 

LISWET8 

100 

106 

2.5747454E-01 

0 

94 

417 

DP 

0 

417 

LISWET9 

100 

103 

2.1543892E+01 

0 

28 

263 

DP 

0 

263 

LISWETIO 

100 

106 

2.5874831E-01 

0 

68 

378 

DP 

0 

378 

LISWETll 

100 

106 

2.5704145E-01 

0 

68 

379 

DP 

0 

379 

LISWET12 

100 

106 

9.1994948E+00 

0 

37 

460 

DP 

0 

460 

LOTSCHD 

7 

12 

2.3984158E+03 

4 

8 

16 

DP 

0 

16 

MOSARQPl 

10 

100 

-1.5420010E+02 

0 

102 

52 

DP 

0 

52 

M0SARqP2 

10 

100 

-2.0651670E+02 

0 

100 

33 

DP 

0 

33 

NASH 

24 

72 

infeasible 

5 

5* 

24* 

DP 

0 

24 

OBSTCLAE 

1 

529 

1.6780270E+00 

0 

605 

178 

DP 

0 

178 

DBSTCLAL 

1 

529 

1.6780270E+00 

0 

263 

263 

PD 

263 

0 

QBSTCLBL 

1 

529 

6.5193252E+00 

0 

469 

469 

PD 

469 

0 

QBSTCLBM 

1 

529 

6.5193252E+00 

0 

484 

189 

DP 

0 

189 

QBSTCLBU 

1 

529 

6.5193252E+00 

0 

303 

303 

PD 

303 

0 

DSLBqP 

1 

8 

6.2500000E+00 

0 

6 

0 

PD 

0 

0 

PENTDI 

1 

500 

-7.5000000E-01 

0 

2 

2 

PD 

2 

0 

P0WELL20 

100 

100 

5.2703125E+04 

49 

52 

99 

DP 

0 

99 
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Table 1: Results for PDQP and SQOPT on 121 CUTEst QPs. (continued) 


Name 

m 

n 

Objective 

SQOPT 

Phsl Itn 

Itn 

PDQP 

Order P-Itn 

D-Itn 

PRIMAL1 

85 

325 

-3.5012967E-02 

0 

217 

70 

PD 

70 

0 

PRIMAL2 

96 

649 

-3.3733671E-02 

0 

407 

97 

PD 

97 

0 

PRIMAL3 

111 

745 

-1.3575583E-01 

0 

1223 

102 

PD 

102 

0 

PRIMAL4 

75 

1489 

-7.4609064E-01 

0 

1264 

63 

PD 

63 

0 

PRIMALCl 

9 

230 

-6.1552516E+03 

0 

18 

5 

PD 

5 

0 

PRIMALC2 

7 

231 

-3.5513063E+03 

0 

3 

5 

PD 

5 

0 

PRIMALC5 

8 

287 

-4.2723256E+02 

0 

10 

6 

PD 

6 

0 

PRIMALC8 

8 

520 

-1.8309432E+04 

0 

30 

6 

PD 

6 

0 

QPCBLEND 

74 

83 

-7.8425425E-03 

0 

111 

182 

PD 

182 

0 

qPCBOEIl 

351 

384 

1.1503952E+07 

415 

1055 

793 

PD 

395 

398 

qPCB0EI2 

166 

143 

8.1719635E+06 

142 

315 

340 

PD 

163 

177 

qPCSTAIR 

356 

467 

6.2043917E+06 

210 

433 

970 

PD 

645 

325 

qUDLIN 

1 

420 

-8.8290000E+06 

0 

419 

419 

PD 

419 

0 

RDW2D51F 

225 

578 

1.1209939E-03 

29 

29 

217 

DP 

0 

217 

RDW2D51U 

225 

578 

8.3930032E-04 

14 

16-f 

219 

DP 

0 

219 

RDW2D52B 

225 

578 

1.0947648E-02 

349 

488 

316^ 

DP 

0 

314 

RDW2D52F 

225 

578 

1.0491239E-02 

29 

191 

414/ 

DP 

0 

414 

RDW2D52U 

225 

578 

1.0455316E-02 

15 

318-f 

219 

DP 

0 

219 

3268 

5 

5 

7.2759576E-12 

0 

8 

0 

PD 

0 

0 

SIM2BqP 

1 

2 

O.OOOOOOOE+00 

0 

1 

1 

PD 

1 

0 

SIMBqP 

1 

2 

6.0185310E-31 

0 

2 

1 

PD 

1 

0 

STCqPl 

30 

65 

4.9452085E+02 

8 

53 

20 

DP 

0 

20 

STCqP2 

128 

257 

1.4294017E+03 

80 

215 

73 

DP 

0 

73 

STEENBRA 

108 

432 

1.6957674E+04 

14 

89 

177 

PD 

2 

175 

TAME 

1 

2 

3.0814879E-33 

0 

1 

1 

PD 

1 

0 

TORSIONl 

1 

484 

-4.5608771E-01 

0 

256 

256 

PD 

256 

0 

T0RSI0N2 

1 

484 

-4.5608771E-01 

0 

544 

144 

DP 

0 

144 

TORSIONS 

1 

484 

-1.2422498E+00 

0 

112 

112 

PD 

112 

0 

T0RSI0N4 

1 

484 

-1.2422498E+00 

0 

689 

288 

DP 

0 

288 

TORSIONS 

1 

484 

-2.8847068E+00 

0 

40 

40 

PD 

40 

0 

TORSIONS 

1 

484 

-2.8847068E+00 

0 

708 

360 

DP 

0 

360 

TORSIONA 

1 

484 

-4.1611287E-01 

0 

272 

272 

PD 

272 

0 

TORSIONB 

1 

484 

-4.1611287E-01 

0 

529 

128 

DP 

0 

128 

TORSIONC 

1 

484 

-1.1994864E+00 

0 

120 

120 

PD 

120 

0 

TORSIDND 

1 

484 

-1.1994864E+00 

0 

681 

280 

DP 

0 

280 

TORSIONE 

1 

484 

-2.8405962E+00 

0 

40 

40 

PD 

40 

0 

TORSIONF 

1 

484 

-2.8405962E+00 

0 

761 

360 

DP 

0 

360 

UBHl 

60 

99 

1.1473520E+00 

11 

40 

112 

DP 

0 

112 

YAO 

20 

22 

2.3988296E+00 

0 

2 

20 

DP 

0 

20 

ZECEVIC2 

2 

2 

-4.1250000E+00 

0 

4 

5 

PD 

5 

0 




i = infeasible, f = 

failed 





8. Summary and Conclusions 

A pair of two-phase active-set methods, one primal and one dual, are proposed for 
convex quadratic programming. The methods are derived in terms of a general 
framework for solving a convex quadratic program with general equality constraints 
and simple lower bounds on the variables. In each of the methods, the search di¬ 
rections satisfy a KKT system of equations formed from Hessian and constraint 
components associated with an appropriate column basis. The composition of the 
basis is specified by an active-set strategy that guarantees the nonsingularity of 
each set of KKT equations. In addition, a combined primal-dual active set method 
is proposed in which a shifted dual QP is solved for a feasible point for the primal 
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Figure 1: Performance profile of number of 

iterations for PDQP and SQOPT on 121 CUTEst 
QP problems. 



Figure 2: Outperforming factors for total 

iterations for each of the 121 CUTEst QP 
problems solved using PDQP and SQOPT. 


(or vice versa), thereby avoiding the need for an initial feasibility phase that ignores 
the properties of the objective function. This approach provides an effective method 
for finding a dual-feasible point when the QP is convex but not strictly convex. Pre¬ 
liminary numerical experiments indicate that this combined primal-dual active set 
method can be an efficient and reliable alternative to conventional two-phase active- 
set methods. Future work will focus on the application of the proposed methods 
to situations in which a series of related QPs must be solved, for example, in se¬ 
quential quadratic programming methods and methods for mixed-integer nonlinear 
programming. 
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A. Appendix 

The appendix concerns some basic results used in previous sections. The first result 
shows that the nonsingularity of a KKT matrix may be established by checking that 
the two row blocks (H ) and (A — M) have full row rank. 


Proposition A.l. Assume that H and M are symmetric, positive semidefinite ma¬ 
trices. The vectors u and v satisfy 


if and only if 



(A.l) 


(A.2) 


Proof. If (A.2) holds, then (A.l) holds, which establishes the “if” direction. Now 
assume that u and v are vectors such that (A.l) holds. Then, 


u'^Hu — u^AA = 0, and v'^Au + v'^Mv = 0. 
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Adding these equations gives the identity u^Hu + v^Mv = 0. But then, the sym¬ 
metry and semidefiniteness of H and M imply vf^Hu = 0 and v^Mv = 0. This can 
hold only if Hu = 0 and Mv = 0. If Hu = 0 and Mv = 0, (A.l) gives A^v = 0 and 
Au = 0, which implies that (A. 2 ) holds, which completes the proof. | 


The next result shows that when checking a subset of the columns of a symmetric 
positive semidefinite matrix for linear dependence, it is only the diagonal block that 
is of importance. The off-diagonal block may be ignored. 


a symmetric, positive semidefinite matrix partitioned 

/Fn H,A 

Us ^ 227 ' 

= if and only if Hnu = 0 . 

Proof. If H is positive semidefinite, then Hu is positive semidefinite, and it holds 
that 

\uy \H12J \^i 2 ^22/ \^J 

if and only if 

“ = £)(o) 

if and only if Huu = 0 , as required. | 


Proposition A.2. Let H be 


as 


Then, 


(S) “' 


In the following propositions, the distinct integers k and I, together with integers 
from the index sets B and J\f define a partition of X = {1, 2, ..., n}, i.e., X = 
BU {fe} U {/} UAf. If w is any n-vector, the ns-vector Wb and tcjv-vector Wn denote 
the vectors of components of w associated with B and N. For the symmetric Hessian 
H, the matrices Hbb and FIjviv denote the subset of rows and columns of H associated 
with the sets B and M respectively. The unsymmetric matrix of components hij 
with i G B and j G M will be denoted by Hbn- Similarly, Ag and Ajv denote the 
matrices of columns associated with B and N. 

The next result concerns the row rank of the (A — M) block of the KKT 
matrix. 


Proposition A.3. If the matrix (ai Uk Ag — M) has full row rank, and there 
exist Axi, Axk, Axg, and Ay such that aiAxi + a^Axh + AgAxg + MAy = 0 with 
Axk / 0, then (^ai Ag — M) has full row rank. 

Proof. It must be established that u"^ {ai Ag —M) = 0 implies that u = 0. For 
a given u, let 7 = —u'^Ok, so that 

{ u ^ 7) = (0 0 0 0). 
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Then, 


0 = [ u ^ 7) ^ 

As Axk / 0, it must hold that 7 = 0, in which case 



( I^xi\ 

(ai Ok Ag -M\ 

Axk 

\ 1 ) 

AXg 


\ -^y J 


1 Axk- 


u^{ ai Uk Ag - M ) = 0. 

As (ai Ofc Ag — M) has full row rank by assumption, it follows that u = 0 and 
(ai Ag — M) must have full row rank. | 


An analogous result holds concerning the (H 


A ^) block of the KKT matrix. 


Proposition A. 4. If ( Hgg A^ ) has full row rank, and there exist quantities Axi^, 
Axg, Ay, and Az^ such that 


(hi. 

hlk al 1 \ 

AXg 

\hgN 

Hgg Al ) 

-Ay 



\-Azk) 


with Azk 7 ^ 0 , then the matrix 

(Kk 

\^Bk ^BB -^bJ 



(A.3) 


has full row rank. 


Proof. Let (/i ) be any vector such that 



As Azk 7 ^ 0 by assumption, it must hold that /i = 0. The full row rank of ( Hgg ) 
then gives v = 0 and 


h Tj\T 


hlk 

Hgg Al 


must have full row rank. Proposition A.l implies that this is equivalent to 


(Kk hlk 

\hgk Hgg Agj 
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having full row rank. | 

The next proposition concerns the primal subiterations when the constraint in¬ 
dex k is moved from B to N. In particular, it is shown that the Ki matrix is 
nonsingular after a subiteration. 

Proposition A.5. Assume that {Axi, Axk, Axb, —Ay, —Azi) is the unique solu¬ 
tion of the equations 


kl 


h 

ai 

VI 

and that Ax^ 7 ^ 0. Then, the matrices Ki and are nonsingular, where 
'h,. 


Kl hli af 


( AxK 


/o\ 

^kk ^Bk 


Axk 


0 

^Bk ^BB 


Axb 

= 

0 

ak Ab -M 


-Ay 


0 

-V 


V / 


Vv 


(A.4) 



Kl 

K \ 

/^kk 

Kk 

K' 

Hbb 


and Kk = h^,. 

Hbb 

Al 

Ab 

-mJ 

V 

Ab 

-M, 


Proof. By assumption, the equations (A.4) have a unique solution with Axk 7 ^ 0. 
This implies that there is no solution of the overdetermined equations 



^kl 

Kl 

T 

^ki 

^kk 

Kk 

K 

^Bl 

^Bk 

Hbb 

Al 

ai 

a-k 


-M 

1 

v 

1 




1\ 



f Axi'^ 
Axk 
Axb 
-Ay 
\-Azi / 


/0\ 

0 

0 

0 

1 

VO/ 


(A.5) 


Given an arbitrary matrix D and nonzero vector /, the fundamental theorem of 
linear algebra implies that if Dw = f has no solution, then there exists a vector v 
such that 7 ^ 0. The application of this result to (A.5) implies the existence of a 


nontrivial vector {Axi, Axk, Axb, —Ay, —Azi, —Azk) such that 


^ K 

^kl 

Kl 

T 

K 

1 \ 

^ki 

^kk 

Kk 

K 

1 

^Bl 

^Bk 

Hbb 

Al 


ai 

VI 

ak 

Ab 

-M 

-1 / 


/ Axi\ 

Axk 

Axb 

-Ay 

-Azi 

\-AzkJ 


/ 0 \ 

0 

0 

0 

VO/ 


(A. 6 ) 


with Azi 0. The last equation of (A. 6 ) gives Axi -|- Azi = 0, in which case 
AxiAzi = —Azf < 0 because Azi 7 ^ 0. Any solution of (A. 6 ) may be viewed as 
a solution of the equations HAx — A^Ay — AT = 0 , AAx + MAy = 0 , Azb = 0 , 
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and Axi = 0 for i G {1, 2, ..., n}\ {/} \ {k}. An argument similar to that used to 
establish Proposition |2.2| gives 

AxiAzi + AxkAzk > 0, 


which implies that Ax^Az^ > 0, with Ax^ 7 ^ 0 and Az^ ^ 0. 

As the search direction is uniqu e, it f ollows from (A.4) that ( H^,. Hgg ) 
has full row rank, and Proposition A.2 implies that (Hgg Ag^ has full row rank. 
Hence, as Azi ^ 0, it follows from (A.6) and Proposition |A.4| that the matrix 


h„ h 


kl 


hh 


^Bl ^Bk ^BB -^B 

has full row rank, which is equivalent to the matrix 

'h. 


h Tjj H nn 


Al 


having full row rank by Proposition A.2 


Again, the search direction is unique and (A.4) implies that (a^ Ag — M) 

implies that (ai As — M) must 


has full row rank. As Axk ^ 0, Proposition 


A.3 


have full row rank. Consequently, Proposition A.l implies that Ki is nonsingular. 

As Axk, Axi, Azk and Azi are all nonzero, the roles of k and I may be reversed 
to give the result that Kk is nonsingular. | 


The next proposition concerns the situation when a constraint index k is moved 
from AA to H in a dual subiteration. In particular, it is shown that the resulting 
matrix Kg defined after a subiteration is nonsingular. 


Proposition A.6. Assume that there is a unique solution of the equations 



^kl 

hsi 

aj 1 

\ 


/ 

Axi\ 


/o\ 

^kl 

^kk 

hi 

al 

1 


Axk 


0 

^Bl 

^Bk 

Hbb 

A'^ 




Axb 


0 

ai 


Ab 

-M 




-Ay 


0 

1 



-1 




-Azi 


1 

V 

1 



/ 


V 

-Azk) 


VO/ 


with Azk A 0- Then, the matrices Ki and Kk are nonsingular, where 


(hii 

h^ 

"'Bl 


. fhkk 

hi 


II 

Hbb 

Al 

1 , and Kk = hg,^ 

Hbb 

Al 

\ai 

Ab 

-m) 

\ 

Ab 

-M, 


(A.7) 


Proof. As (A.7) has a unique solution with Azk A 0) there is no solution of 

1\ 


^hii 

hki 

hi 

T 

hki 

hkk 

hi 

al 

hBi 

hBk 

Hbb 

Al 

ai 

ak 

Ab 

-M 


-1 


/ Axi\ 

Axk 
Axb 
—Ay 

V J 


/o\ 

0 

0 

0 

1 

Vo/ 


(A.8) 
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The fundamental theorem of linear algebra applied to (A. 8 ) implies the existence of 
a solution of 



^kl 

hi 

T 

1 \ 

^kl 

^kk 

hlk 


1 


^Bk 

Hbb 

Al 


ai 

VI 


Ab 

-M 

-1 / 


/ 

Axk 

Axb 

-Ay 

-Azi 

\-AzkJ 


/0\ 

0 

0 

0 

voy 


(A.9) 


with Azi 7 ^ 0. It follows from (A.9) that Axi + Azi = 0. As Azi 7 ^ 0, this 


implies AxiAzi < 0. The solution of (A.9) may be regarded as a solution of the 
homogeneous equations HAx — A^Ay — Az = 0, AAx + MAy = 0, with Azi = 0, 
for i £ B, and Axi = 0, for i G {1,..., n} \ {k} \ {/}. Hence, Proposition 2.2 gives 

AxiAzi + AxkAzk > 0 , 

so that AxkAzk > 0. Hence, it must hold that Axk 7^ 0 and Azk 7^ 0. 

As Axk 7^ 0, Axi 7^ 0, Azk 7^ 0 and Azi 7^ 0, the remainder of the proof is 


analogous to that of Proposition A.5 


The next result gives expressions for the primal and dual objective functions in 
terms of the computed search directions. 


Proposition A.7. Assume that (x, y, z) satisfies the primal and dual equality con¬ 
straints 

Hx + c — A^y — z = 0, and Ax + My — 6 = 0. 


Consider the partition {1, 2, n} = B U {/} U J\f such that Xn + Qn = 0 and 
Zb + rg = 0. If the components of the direction {Ax, Ay, Az) satisfy (2.8), then 
the primal and dual objective functions for (PQP^ ,.) and (DQP^ i.e.. 


fp{x,y) = ^x Hx^y Myc xr x 
fn{x, y, z) = -\x'^Hx - ly^My + b^y - q^z, 


satisfy the identities 

fp{x -\- aAx,y -\- aAy) = fp{x,y) -\- Axfizi + n)a + \ AxiAzia^, 
fnix + aAx,y + aAy,z + aAz) = fD{x,y,z) - Azfixi + qi)a - \AxiAzia^. 


Proof. The directional derivative of the primal objective function is given by 
{Ax'^ Ay'^)Vfp{x,y) = {Ax'^ ^ '^) 

= {Ax^ 

= {AAx + MAy)'^y + Ax'^{z + r) = Axfizi + r;), (A.10b) 
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where the identity Rx + c = AJ'y + 2 ; has been used in (A.10a) and the identities 
AAx + MAy = 0, Axn = 0 and Zg + = 0 have been used in (A. 10b). 

The curvature in the direction (Z\x, Ay) is given by 


Ay'^)V‘^fp{x,y) ^v'^) 


fAx\ 
MJ [AyJ 


= AxiAzi, 


(ATI) 


where the last equality follows from Proposition |2.2[ 

The directional derivative of the dual objective function is given by 

/ -Hx \ 

i^Ax'^ Ay"^ Az^) Vf d{x, y, z) = [Ax’^ Ay"’" Az’") j —My + 6 1 (A.12a) 

=—Ax’^Hx + Ay'’"{—My+ b) — Az’" q (A.12b) 

= —{A^Ay + Az)"’'x + Ay"’"{—My + b) — Az’" q 

(A.12c) 

= —Ay"'"{Ax + My — b) — Az"’"{x + q) (A.12d) 
=-Azi{xi +qi), (A.12e) 


where the identity HAx—A’"Ay—Az = 0 has been used in (A.12c) and the identities 
Ax + My — 6 = 0, Xjv + Q'iv = 0 and Azp = 0 have been used in (A.12e). 

As 2 ; only appears linearly in the dual objective function, it follows from the 
structure of the Hessian matrices of fp{x,y) and fp,{x,y,z) in combination with 
(A.ll) that 

{Ax’" Ay’' Az’)V^fDix,y,z)(Ayj=-{Ax’ Ay’)V‘^fp{x,y)(^^ 


= -AxiAzi. 


The final result shows that there is no loss of generality in assuming that {AM) 
has full row rank in (PQPg ,.). 

Proposition A.8. There is no loss of generality in assuming that {A M) has full 
row rank in (PQPg ,.). 


Proof. Let {x, y, z) be any vector satisfying ( |2.1a[ ) and (2.1b). Assume that 
{A M) has linearly dependent rows, and that {A M) and 6 may be partitioned 
conformally such that 


(.4 A/) = ( 


Ai Mil Mi2 

A2 ^12 ^22 


and 6 = 


with (Ai Mil Mi 2 ) having full row rank, and 

(T2 Mfa M22) = N {Ai Mil M12) , 


(A.13) 
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with Ai G and A 2 G for some matrix N G ]R™-2xmi^ From the linear 

dependence of the rows of M), it follows that x, y and z satisfy (2.1a) and 
(2.1b) if and only if 


Hx + c - A^y^ - Aly 2 - z = 

Aix + Miiyi + Mi 2 y 2 — 61 = 0 and 62 = Nbi. 


It follows from (A.13) that M 12 = MniV’^ and = AJn'^, so that x, y and z 
satisfy (2.1a) and (2.1b) if and only if 


Hx + c- AJ{y^ + - z = 0, 

Aix + Mii(yi + N'^y 2 ) — 61 = 0 and 62 = Nbi. 


We may now define yi 


yi + A^''y2 and replace (2.1b) and (2.1a) by 




Hx + c — A{yi — z = 0, 
Aix + Miiyi - 61 = 0. 


A.2 


implies that 


By assumption, (Ai Mn M12) has full row rank. Proposition 
( Ai Mil ) has full row rank. This gives an equivalent problem for which (Ai Mu ) 
has full row rank. I 













